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COLISEUM  is  a  application  framework  that  integrates  plasma  propagation  schemes  and 
arbitrary  3D  surface  geometries.  Using  Particle-in-Cell  (PIC)  schemes  to  model  the  plasma 
propagation  high  fidelity  modeling  of  the  plasma  and  its  interaction  with  the  surfaces  is 
possible.  In  order  to  improve  the  computational  performance  of  the  Particle-in-Cell  (PIC) 
scheme  within  COLISEUM  (AQUILA),  accelerate  techniques  have  been  developed  that 
significantly  decrease  the  amount  of  CPU  time  needed  to  obtain  a  steady-state  solution. 
These  schemes  have  been  demonstrated  to  decrease  the  CPU  time  from  3  to  24  times 
with  little  appreciable  differences  in  the  global  particle  properties  and  number  densities. 
This  works  investigates  the  differences  in  the  local  plasma  properties  that  result  from  the 
application  of  the  different  acceleration  techniques.  In  particular,  the  number  densities  and 
velocity  distributions  of  the  ions  and  neutrals  demonstrate  that  the  solution  acceleration 
schemes  produce  very  similar  solutions  outside  of  the  main  path  of  the  plasma  source. 
Within  the  main  path  of  the  plasma  source  the  local  plasma  properties  show  marked 
differences  that  might  be  associated  with  the  time  steps  associated  with  these  schemes 
and/or  the  collision  modeling  scheme  within  AQUILA. 


Nomenclature 


c  Velocity,  [m/s] 

Ci  Three-dimensional  velocity  space,  [m/s  x  m/s  x  m/s] 
cr  Relative  velocity  between  two  particles  [m/s] 

/  Velocity  distribution  function,  [-] 

Mass  flow  rate,  [kg/s] 

Number  of  computational  particles 
n  Particle  number  density  [m~3] 

P  Probability  of  collision 

S  Surface  of  integration 

t  Time,  [s] 

Te  Electron  Temperature  [eV] 

V  Volume  of  computational  cell  [m3] 

Wp  Ratio  of  physical  particles  to  computational  particles 

x  Position,  [m] 

Subscripts 
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max  Maximum  value 

s  Surface  property 

Conventions 

NTC  No- Time-Counter  Collision  Model 
SCN  COLISEUM  cases  run  without  subcyling 
SCY  COLISEUM  cases  run  with  subcyling 
VDF  Velocity  Distribution  Function 

Symbols 

(j)  Electrostatic  potential  [V] 

a t  Total  collision  cross-section  [m2] 

r  Characteristic  time  [s] 


I.  Introduction 

The  Air  Force  Research  Laboratory  has  developed  an  application  framework  that  integrates  plasma  prop¬ 
agation  schemes  with  arbitrary  3D  surface  geometries1  in  order  to  investigate  the  interactions  between 
the  plasma  plume  from  electric  propulsion  thrusters  and  the  spacecraft.  A  hybrid  Particle-in-Cell  plasma 
model  within  COLISEUM,  called  AQUILA,2  is  the  basis  of  this  work.  The  COLISEUM  framework  is  used  to 
model  the  thrusters  in  space  environment  for  prediction  of  the  interactions  and  also  to  model  the  thrusters 
in  vacuum  chambers  in  order  to  validate  the  models  used  in  COLISEUM  against  experimental  data.  In 
order  to  improve  the  computational  performance  of  the  AQUILA  model,  acceleration  techniques  have  been 
developed  that  significantly  decrease  the  amount  of  CPU  time  needed  to  get  the  simulation  to  a  steady-state. 
These  schemes  have  been  demonstrated  to  decrease  the  CPU  time  by  up  to  24  times  with  little  appreciable 
differences  in  the  global  particle  properties. 

Two  previous  studies  have  been  performed  on  the  acceleration  schemes  within  COLISEUM.  Gibbons  et 
al.3  also  demonstrated  the  subcycling  technique  in  which  the  slower  moving  neutrals  are  propagated  at  a 
larger  time  step  than  the  faster  moving  ions.  Gibbons  et  al.4  demonstrated  the  use  of  the  subcycling  scheme 
with  a  scheme  that  decoupled  the  modeling  of  surface  interactions  from  the  plume  propagation.  It  utilizes 
particle  sources  from  the  surfaces  in  place  of  the  self-consistent  surface  interaction  modeling.  Both  of  these 
schemes  are  intended  to  provide  increased  convergence  rates  to  a  steady-state  solution.  The  desire  is  to  end 
up  with  the  same  final  plasma  distribution  with  or  without  the  acceleration  techniques.  The  present  study 
intends  to  provide  a  detailed  analysis  of  the  resulting  plasma  properties  and  any  differences  in  the  local 
plasma  properties  that  result  from  the  application  of  the  different  acceleration  techniques. 

II.  Computational  Techniques 

This  section  provides  a  brief  overview  of  the  major  computational  techniques  within  COLISEUM.  Consult 
the  User’s  Manual  for  more  details.5,6 

A.  COLISEUM  Framework 

COLISEUM  is  a  computational  framework  that  provides  a  tool  that  can  be  used  to  model  the  interaction 
of  plasmas  with  arbitrary  surfaces  in  three-dimensional  space.3  These  simulations  can  be  of  a  plasma  in  a 
contained  domain  or  in  an  open  domain.  This  allows  the  simulation  of  experiments  in  vacuum  chambers  as 
well  as  plasmas  in  the  low  density  space  environment.  The  primary  focus  of  COLISEUM  is  to  investigate 
the  erosion  associated  with  the  plasma  particles  impacting  surfaces,  known  as  sputtering.  This  sputtered 
material  may  be  re-deposited  on  other  surfaces.  COLISEUM  provides  the  integration  of  the  CAD  surface 
modeling,  the  plasma  propagation,  and  the  sputtering  of  material  within  one  consistent  framework. 

1.  Surface  Modeling 

The  surface  model  can  be  input  in  a  variety  of  standard  formats  including  ANSYS  and  ABAQUS  formats. 
In  addition,  surface  properties  are  also  specified  in  order  to  differentiate  the  various  materials  that  may 
compose  each  surface.  In  order  to  accurately  model  the  sputtering,  additional  material  information  must  be 
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specified  that  describes  the  interaction  between  particles  that  may  be  impacting  surfaces  and  the  material 
that  the  surfaces  are  composed  of.  This  is  known  as  the  material  interaction  parameters  in  COLISEUM. 


2.  Plasma  Modeling 


The  plasma  modeling  within  COLISEUM  has  two  major  components  to  it.  The  first  is  the  modeling  of  the 
source.  Sources  are  surfaces  within  the  geometry  that  particles  will  be  emitted.  To  model  sources,  a  velocity 
distribution  function,  VDF,  must  be  specified  throughout  the  surface  of  the  source.  In  general,  the  VDF  is 
a  function  of  space,  time,  and  obviously  velocity,  /  (x,  c,  t).  This  is  used  to  determine  the  mass  flow  rate,  m, 
as 


rhs 


S  J  ci 


fs  (x,  c,  t)  dCidS 


(1) 


where  Cj  is  the  three-dimensional  velocity  space,  and  S  is  the  surface  of  the  source.  The  VDF  does  not 
need  to  specified  direction,  as  COLISEUM  provides  surface  modules  that  model  common  sources:  (1)  user 
specified  flux  information  typically  from  experimental  data,  (2)  user  specified  flux  and  velocity  information 
typically  from  experimental  data,  and  (3)  a  shifted  Maxwellian  distribution  with  the  velocity  shift  normal 
to  the  surface. 

The  second  major  component  of  the  plasma  modeling  within  COLISEUM  is  the  plasma  simulation  itself. 
COLISEUM  was  developed  with  the  idea  of  supporting  any  number  of  plasma  modeling  schemes.  One  of 
the  original  models  is  the  prescribed  plume  model  which  simply  imports  the  previously  obtained  plasma 
properties  of  the  flow  field.  The  particle  fluxes  on  the  surfaces  are  then  directly  obtained  from  this  data. 
The  other  original  model  is  the  ray  tracing  model  which  traces  the  particle  trajectory  from  the  sources 
without  accounting  for  the  electrostatic  potential  field  forces.  The  particle  fluxes  on  the  surfaces  are  then 
able  to  be  determined. 

As  COLISEUM  has  matured,  more  sophisticated  plasma  modeling  modules  have  been  developed.  DRACO 
from  Virginia  Tech7  is  a  Cartesian  cell  based,  finite-element  PIC-DSMC  simulation.  AQUILA  from  MIT8 
is  an  unstructured  tetrahedral  cell  based,  finite-element  PIC-DSMC  simulation.  It  is  AQUILA  that  is  being 
used  as  the  basis  of  this  investigation. 


3.  Sputtering  Modeling 

The  sputtering  models  in  COLISEUM  are  based  on  standard  models  from  Roussel  et  al.9  and  Gardner  et 
ah, 10  Kannenberg  et  al.,11  and  Yamamura  et  al.12  Coupling  the  sputtering  models  with  the  re-deposition 
process  has  been  included  in  order  to  account  for  how  the  re-deposited  material  may  itself  induce  sputtering. 
This  allows  a  more  accurate  model  of  the  final  surface  deposition  characteristics. 


B.  Particle  Propagation  Scheme 

The  time  integration  scheme  used  within  AQUILA  to  propagate  the  plasma  particles  is  the  standard  leap  frog 
scheme13  which  is  second  order  accurate  in  time.  The  electrostatic  forces  are  modeled  using  the  electrostatic 
potential  equation  with  the  inclusion  of  space-charge  effects.8  A  finite  element  formulation  is  used  to  solve 
the  potential  equation  with  a  Newton-method  type  scheme  to  handle  the  nonlinear  nature  of  the  resulting 
equations. 


C.  Collision  Modeling 

The  collision  modeling  within  AQUILA14  is  based  on  the  No-Time-Counter,  NTC,  method  of  Bird.15  The 
probability  of  a  collision  between  two  particles  is  given  as 

=  Wp  ( CTTCr )  At 

V  1  ’ 

where  Wp  is  the  ratio  of  physical  particles  to  computational  particles,  aT  is  the  total  collision  cross-section, 
cr  is  the  relative  speed  between  the  two  particles,  and  V  is  the  volume  of  the  computational  cell  containing 
particles.  Similarly,  the  maximum  probability  of  a  collision  is 


P  = 

1  m  py 


WP  Qrcr)max  At 
V 


(3) 
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The  NTC  scheme  samples  only  a  fraction  of  the  total  number  of  particle  pairs  in  the  computational  cell, 
and  similarly  increases  the  probability  of  collision  of  each  sampled  particle  pairs.  Within  COLISEUM,  only 
Pm&xNpNq  particle  pairs  are  chosen  from  species  p  and  q.  Thus,  the  lower  the  maximum  probability  of 
collision,  the  fewer  collision  samples  are  taken.  The  resulting  collision  probability  for  a  sampled  collision 
pair  is  then 


P  = 


(It  cr 

(aTCr)max 


(4) 


A  simple  accept/reject  scheme  can  be  used  on  this  probability.  Notice  that  this  scheme  will  sample  the 
appropriate  number  of  collision  pairs  only  when  an  accurate  maximum  probability  has  been  determined. 
Therefore,  this  scheme  will  produce  accurate  collision  rates  only  after  a  large  number  of  collisions  have 
occurred  so  that  the  maximum  probability  term  has  been  reasonably  determined. 


D.  Subcycling  Scheme 

The  subcycling  scheme  within  COLISEUM  utilizes  the  fact  that  there  is  a  significant  difference  between  the 
collision  times  scales  and  the  ion  characteristic  times.  The  ion  characteristic  time  is  based  on  the  spatial 
resolution  of  the  computational  grid  which  is  being  used  to  model  the  electrostatic  potential  field.  In  order 
to  keep  the  electrostatic  forces  on  the  particles  varying  smoothly  as  they  travel  through  the  domain,  the 
ions  should  not  travel  more  than  a  third  of  a  cell  within  one  time  step.  For  the  simulation  to  be  modeled 
in  this  paper,  the  ion  velocity  is  20  km/s  and  the  neutral  velocity  is  200  m/s.  For  characteristic  length 
of  the  smallest  volume  as  0.01  m  this  results  in  the  ion  characteristic  time  of  5  x  10-7  s  and  the  neutral 
characteristic  time  of  5  x  10-5  s.  This  is  a  factor  of  100  difference  between  the  two. 

Using  the  following  relations  for  the  elastic  collisions  between  neutrals  and  neutral  ions16  to  characterize 
the  collision  cross-sections 


_el  _  2.117  x  10~18 

aXe-Xe  ~  o.24 

el  8.2807  x  10"16 

aXe-Xe+  ~  “ 

and  for  the  charge-exchange  collisions  between  neutrals  and  ions17 

aXe-Xe+  =  L1872  x  10“2°  [-23.3  log  (Cr)  +  188.81] 


(5) 

(6) 

(7) 


a  mean  time  to  collision  can  be  determined  as 


1 

r  =  - 

ncrcr 


(8) 


Table  1  shows  the  characteristic  times  for  the  simulation  to  be  modeled  in  this  paper  and  a  maximum  particle 
number  density  of  1018  m'3.  Included  in  these  calculations  is  the  high  speed  neutrals  that  will  result  from 
previous  charge  exchange  collisions.  Thus,  the  neutrals  in  the  collision  calculations  could  have  the  low  speed 
200  m/s  value  or  the  high  speed  20000  m/s  value.  Clearly  a  simulation  time  step  less  than  4.7  x  10-5  s  is 
needed  to  resolve  the  collision  time  scales. 


Table  1.  Collision  Characteristic  Times 


Collision  Type 

cr  [m/s] 

a  [m2] 

T  [S] 

Xe-Xe  Elastic 

200 

5.94  x  10"19 

8.42  x  10"3 

Xe-Xe  Elastic 

10000 

2.32  x  10~19 

4.28  x  10"4 

Xe-Xe  Elastic 

20000 

1.97  x  10"19 

2.54  x  10"4 

Xe-Xe+  Elastic 

10000 

8.20  x  10"2° 

1.21  x  10"3 

Xe-Xe+  Elastic 

20000 

4.14  x  10"2° 

1.21  x  10"3 

Xe-Xe+  Charge  Exchange 

10000 

1.13  x  10~18 

8.73  x  10"5 

Xe-Xe+  Charge  Exchange 

20000 

1.05  x  10"18 

4.75  x  10"5 
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Therefore,  there  is  only  one  physical  phenomena  that  requires  a  time  step  in  the  order  of  10  1  s,  and 


that  is  capturing  the  electrostatic  forces  applied  to  the  ions. 

Computing  one  complete  computational  cycle  encompasses  the  following  steps 

1.  Subcycle  Fast  Particles 

(a)  Move  Fast  Particles 

(b)  Inject  Fast  Particles 

(c)  Update  Electrostatic  Fields 

2.  Move  Slow  Particles 

3.  Inject  Slow  Particles 

4.  Perform  Collisions 

By  only  moving  the  slow  particles  a  fraction  of  the  number  times  that  the  fast  particles  must  be  moved 
as  well  as  performing  the  collisions  on  the  coarse  time  step  a  significant  amount  of  computational  effort  is 
saved. 

E.  Velocity  Distribution  Function  Probe 

In  order  to  determine  the  local  plasma  properties,  a  new  probe  was  introduced  into  the  COLISEUM  probe 
architecture.  This  probe  samples  a  region  is  space  (currently  a  computational  cell)  and  stores  the  velocity 
of  every  particle  of  a  specified  type  that  resides  within  the  cell.  The  frequency  of  sampling  can  be  adjusted 
as  well  as  the  start  and  end  times  of  the  sampling.  Once  the  sampling  is  completed,  the  probe  sorts  the 
particles  into  bins  of  user  specified  sizes.  The  results  can  be  written  out  as  a  table  of  the  non-empty  bins,  or 
as  a  Tecplot  formatted  structured  grid  data  file.  Further  processing  is  possible  with  this  data  if  only  binning 
on  particle  speed  is  desired. 


III.  Results 


The  following  results  are  all  for  the  same  test  problem.  First,  solution  convergence  is  demonstrated  by 
using  a  fine  time  step,  that  is  on  the  order  of  the  characteristic  time  step  of  the  electrostatic  forces,  and 
demonstrating  that  further  refinement  of  the  time  step  does  not  result  in  any  appreciable  change  in  the 
solution.  Second,  a  solution  is  presented  with  a  coarse  time  step,  that  is  on  the  order  of  the  collision  time 
scale,  and  is  compared  to  the  fine  time  step  solution.  Next,  a  solution  is  presented  using  the  subcycling 
scheme  discussed  above,  with  the  ions  moving  at  the  fine  time  step  and  the  neutrals  moving  at  the  coarse 
time  step. 

A.  Test  Problem  Description 

The  test  problem  is  a  highly  simplified  geometry  based  on  a  plasma  source  within  a  vacuum  chamber. 
Figure  1  shows  the  surface  meshes  associated  with  the  test  problem.  The  plasma  source  is  a  small  cylinder 
with  the  cylinder  axis  aligned  with  the  z-axis.  The  plasma  is  emitted  in  the  positive  z-direction  from  that 
particular  face  of  the  cylinder.  The  vacuum  chamber  is  simplified  to  a  cylinder  with  the  cylinder  axis  again 
aligned  with  the  z-axis.  The  plasma  source  is  firing  at  one  end  of  the  cylinder  and  the  other  end  of  the 
cylinder  is  a  particle  sink  such  that  any  particle  that  hits  that  surface  leaves  the  computational  domain. 

The  chamber  has  a  diameter  of  1.5  m  and  a  length  of  2  m.  The  plasma  source  as  a  diameter  of  0.1  m 
and  a  length  of  0.1  m.  The  distance  between  the  plasma  source  face  and  the  chamber  face  is  1.3  m 

The  plasma  source  is  composed  of  two  particle  types,  a  low  speed  neutral  and  a  high  speed  ion.  Both 
are  modeled  using  the  drifting  Maxwellian  source  model  within  COLISEUM.5  The  neutral  drift  velocity  is 
200  m/s  and  temperature  is  700  K.  The  ion  drift  velocity  is  20  km/s  and  temperature  is  10  eV. 

The  electrostatic  potential  is  modeled  using  the  quasi-neutral  model  within  AQUILA.14  This  applies  the 
quasi-neutrality  assumption  and  inverting  Boltzmann’s  equation  to  obtain  an  expression  for  the  electrostatic 


potential 


(9) 
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where  the  electron  temperature,  Te,  is  set  to  2  ev  and  the  reference  electron  number  density  and  potential 
is  specified  to  be  at  a  potential  of  0  V  just  in  front  of  the  thruster  face. 


In  order  to  examine  the  similarities  of  the  local 
properties  of  the  plasma  between  the  three  cases,  the 
velocity  distribution  for  the  neutrals  and  ions  were 
obtained  0.1  m  in  front  of  the  plasma  source  as  well 
as  0.28  m  above  the  thruster  face.  The  first  sampling 
will  examine  the  plume  modeling  capabilities,  while 
the  second  sampling  will  examine  the  capabilities  to 
model  the  plasma  outside  of  the  plume. 

The  computer  that  these  simulations  were  per¬ 
formed  for  the  timing  results  is  a  dual  processor 
AMD  Opteron  242  system  with  2  GB  of  RAM  with 
an  additional  2  GB  of  swap  space.  Use  of  the  ma¬ 
chine  was  minimized  while  the  cases  were  running, 
and  all  cases  resided  in  physical  memory,  so  the  swap 
space  was  not  utilized  except  to  move  other  non- 
essential  applications  out  of  RAM  at  the  start  of 
the  simulation. 

In  order  to  quickly  distinguish  between  the  vari¬ 
ous  cases  to  be  run,  the  acronym  SCN  will  designate 
cases  without  the  use  of  subcycling  and  SCY  will 
designate  the  cases  with  subcycling. 


Figure  1.  Simple  Chamber  Geometry 


B.  Solution  Convergence  Demonstration 

With  the  baseline  fine  time  step  for  this  simulation 

established  as  2.5  x  1CU'  s,  three  cases  were  run  to  demonstrate  the  convergence  of  the  solution  at  this  time 
step.  One  case  was  at  twice  the  baseline  time  step,  5.0  x  10~'  s,  and  one  at  half  the  baseline  time  step, 
1.25  x  10” '  s.  Each  case  was  computed  to  a  final  computational  time  of  0.25  s.  Table  2  shows  the  collision 
rates  for  the  three  cases.  There  is  very  little  difference  between  all  of  the  collision  rates  between  the  three 
cases  with  the  differences  between  the  cases  most  certainly  associated  with  the  statistical  scatter  associated 
with  these  types  of  schemes.  Notice  that  halving  or  doubling  the  step  size  approximately  doubled  or  halved 
the  amount  of  time  required  to  calculate  the  solution. 

Table  2.  Collision  Rates  for  Solution  Convergence  Cases 


Scheme 

Compute 
Time  [hr] 

Time 

Step  [s] 

Total 

Collisions  [#/s] 

Xe-Xe+  Charge 
Exchange  [#/s] 

Xe-Xe 
Elastic  [#/s] 

Xe-Xe+ 
Elastic  [#/s] 

SCN 

107.2 

1.25  x  10"7 

5.301  x  107 

5.046  x  107 

2.519  x  105 

2.293  x  106 

SCN 

58.9 

2.5  x  10-7 

5.325  x  107 

5.069  x  107 

2.560  x  105 

2.303  x  106 

SCN 

30.9 

5.0  x  10-7 

5.310  x  107 

5.056  x  107 

2.541  x  105 

2.287  x  106 

Figure  2  shows  the  evolution  of  the  total  number  of  neutrals  and  ions  as  well  as  the  total  number  of 
particles.  These  counts  are  nearly  identical  with  any  differences  with  the  statistical  scatter. 

Figure  3  shows  the  final  number  density  for  the  neutrals  for  the  three  cases.  Figure  4  shows  the  final 
number  density  for  the  ions  for  the  three  cases,  and  Figure  5  shows  the  final  electrostatic  potential  for  the 
three  cases.  All  of  these  again  show  on  very  minor  differences  between  the  three  cases,  and  these  are  when 
the  number  densities  are  low  and  are  associated  with  the  statistical  scatter.  The  outer  wings  of  the  plume 
are  captured  in  all  three  cases  as  can  be  seen  in  the  ion  number  density  and  the  electrostatic  potential.  The 
time  step  is  sufficient  to  capture  the  ion  curved  trajectories  as  is  shown  behind  the  plasma  source  with  ions 
occupying  some  of  the  region  behind  the  plasma  source.  The  very  low  ion  density  just  behind  the  ion  source 
is  due  to  the  ions  colliding  with  that  surface  and  reflecting  back  as  accommodated  neutrals. 

Next,  the  local  properties  of  the  plasma  will  be  compared  between  the  three  cases.  The  neutral  sampling 
within  in  the  plume  is  shown  in  Figure  6  for  both  the  velocity  distribution  as  well  as  the  speed  distribution, 
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particle  count 


Figure  2.  Global  Particle  Counts  for  Solution  Convergence  Cases 


Figure  3.  Final  Neutral  Number  Density  for  Solution  Convergence  Cases 
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(a)  SCN  At  =  1.25  x  10~7  s 


(b)  SCN  At  =  2.5  x  10~7  s 


(c)  SCN  At  =  5.0  x  10"7  s 


Figure  4.  Final  Ion  Number  Density  for  Solution  Convergence  Cases 
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Figure  5.  Final  Electrostatic  Potential  for  Solution  Convergence  Cases 


where  the  velocity  distribution  is  showing  the  outer  edge  of  the  velocity  space  domain  that  is  populated  with 
particles  and  the  surfaces  are  colored  by  z-velocity  since  that  is  the  primary  velocity  direction.  Figure  7 
shows  the  ion  sampling  within  the  plume.  Once  again,  there  is  very  little  differences  between  the  three 
cases.  A  bimodal  distribution  is  seen  in  the  ion  distribution  with  the  lower  speed  ions  being  the  result  of  the 
charge-exchange  collisions  between  the  neutrals  and  the  ions.  The  corresponding  high  speed  neutrals  can 
be  seen  in  the  neutral  distributions,  but  the  clarity  is  obscured  because  of  the  relative  low  occurrence  of  the 
particles  compared  to  the  computational  particle  weighting  used  for  the  neutrals. 


(a)  SON  Velocity  At  =  1.25  x  10  7  s 


speed  (m/s) 

(d)  SCN  Speed  At  =  1.25  x  10“7  s 


(b)  SCN  Velocity  At  =  2.5  x  10  7  s 


speed  (m/s) 


(e)  SCN  Speed  At  =  2.5  x  10  7  s 


(c)  SCN  Velocity  At  =  5.0  x  10  7  s 


speed  (m/s) 

(f)  SCN  Speed  At  =  5.0  x  lO"7  s 


Figure  6.  Neutral  Velocity  and  Speed  Distributions  0.1  m  in  Front  of  Thruster  for  Solution  Convergence  Cases 


The  neutral  sampling  outside  the  plume  is  shown  in  Figure  8  for  both  the  velocity  distribution  as  well 
as  the  speed  distribution.  Figure  9  shows  the  ion  sampling  outside  the  plume.  In  this  case  there  is  very 
little  difference  between  the  three  cases,  but  some  minor  differences  do  occur.  The  majority  of  the  neutrals 
that  are  being  sampled  in  this  region  are  most  likely  from  the  reflected  from  the  walls.  That  is  why  the 
most  probable  velocity  is  so  low.  However,  the  5.0  x  10~7  s  case  does  show  some  added  noise  in  the  higher 
velocities.  Notice  that  at  this  time  step,  we  are  very  close  to  the  ion  characteristic  time  step,  so  the  increase 
in  the  time  step  size  might  have  altered  the  ion  trajectories  enough  to  alter  the  neutral  distribution  in  the 
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Figure  7.  Ion  Velocity  and  Speed  Distributions  0.1  m  in  Front  of  Thruster  for  Solution  Convergence  Cases 


higher  speeds.  There  still  exists  a  bimodal  distribution  in  the  ion  velocity,  but  the  number  of  low  speed  ions 
is  significantly  less  than  what  were  in  front  of  the  thruster.  Also,  the  most  probable  speed  of  the  ions  has 
dropped  from  the  20  km/s  that  it  was  in  front  of  the  thruster  to  around  4  km/s.  This  is  because  the  ions 
that  make  it  to  this  region  of  the  flow  are  the  low  speed  charge  exchange  ions  that  have  been  accelerated 
through  the  electrostatic  potential  field  into  this  region. 

C.  Coarse  Time  Step  Solution 

Now  the  coarse  time  step  case  can  be  compared  to  the  fine  time  step.  For  this  case  all  particles  and  physical 
processes  are  propagated  at  a  time  step  of  2.5  x  10-5  s.  From  the  previous  discussions,  it  is  known  that  this 
time  step  is  larger  than  the  characteristic  time  associated  with  the  electrostatic  forces,  5.0  x  10~5  s,  and  is 
very  close  to  the  smallest  characteristic  time  associated  with  the  collision  modeling,  4.7  x  10~5  for  the  high 
speed  neutrals  and  ions  charge  exchange  collisions.  Table  3  shows  the  resulting  collision  rates  along  with  the 
fine  time  step  collision  rates. 

Table  3.  Collision  Rates  for  Coarse  Time  Step 

Compute  Time  Total  Xe-Xe+  Charge  Xe-Xe  Xe-Xe + 

Scheme  Time  [hr]  Step  [s]  Collisions  [#/s]  Exchange  [#/s]  Elastic  [#/s]  Elastic  [#/s] 
SCN  58.9  2.5  x  KT7  5.325  x  107  5.069  x  107  2.560  x  105  2.303  x  106 

SCN  1.73  2.5  x  10~5  4.945  x  107  4.726  x  107  2.572  x  105  1.936  x  106 


The  first  thing  to  notice  from  this  figure  is  that  taking  100  times  fewer  time  steps  results  in  a  significant 
decrease  in  compute  time,  by  a  factor  of  around  34.  Unfortunately,  only  the  Xe-Xe  elastic  collision  rate  is 
the  same  between  the  coarse  and  fine  time  steps.  The  elastic  and  charge  exchange  collisions  associated  with 
the  Xe-Xe+  pairs.  Therefore  by  not  adequately  resolving  the  ion  trajectory  there  has  been  a  decrease  in  the 
collision  rate  associated  with  the  ions.  This  could  be  caused  by  the  ions  traveling  entirely  through  the  high 
density  region  in  front  of  the  plasma  source  (where  collisions  are  most  likely)  before  the  ions  can  participate 
in  a  significant  number  of  collision  events.  Taking  the  nominal  ion  velocity  of  20  km/s  and  the  time  step  of 
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(a)  SCN  Velocity  At  =  1.25  x  10  7  s 


speed  (m/s) 

(d)  SCN  Speed  At  =  1.25  x  10"7  s 


(b)  SCN  Velocity  At  =  2.5  x  10  7  s 


speed  (m/s) 

(e)  SCN  Speed  At  =  2.5  x  10— 7  s 


(c)  SCN  Velocity  At  =  5.0  x  10  7  s 


speed  (m/s) 

(f)  SCN  Speed  At  =  5.0  x  10"7  s 


Figure  8.  Neutral  Velocity  and  Speed  Distributions  0.28  m  Above  Thruster  Face  for  Solution  Convergence 
Cases 


(a)  SCN  Velocity  At  =  1.25  x  10  7  s  (b)  SCN  Velocity  At  =  2.5  x  10  7  s 


(d)  SCN  Speed  At  =  1.25  x  10“7  s 


(e)  SCN  Speed  At  =  2.5  x  10“7  s 
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(c)  SCN  Velocity  At  =  5.0  x  10  7  s 


speed  (m/s) 


(f)  SCN  Speed  At  =  5.0  x  10"7  s 


Figure  9.  Ion  Velocity  and  Speed  Distributions  0.28  m  Above  Thruster  Face  for  Solution  Convergence  Cases 
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2.5  x  10-5  s  yields  a  distance  traveled  by  an  ion  of  0.5  m.  This  is  a  significant  distance  from  the  plasma 
source  and  is  certainly  outside  of  the  high  density  region,  compare  this  with  where  the  high  density  regions 
were  in  the  fine  time  step  cases  from  Figures  3  and  4. 

While  significant  differences  are  seen  in  the  collision  rates,  Figure  10  shows  the  evolution  of  the  total 
number  of  neutrals  and  ions  as  well  as  the  total  number  of  particles  is  very  similar.  Thus,  the  differences 
between  the  two  cases  must  be  for  only  a  small,  but  significant,  fraction  of  the  total  particles. 


(a)  SCN  At  =  2.5  x  10"7  s  (b)  SCN  At  =  2.5  x  10~5  s 

Figure  10.  Global  Particle  Counts  for  Coarse  Time  Step 


Figure  11  shows  the  final  number  density  for  the  neutrals  for  the  coarse  and  fine  time  step  cases.  Even 
for  this  case  the  neutral  number  density  is  fairly  similar.  It  appears  that  the  coarse  time  step  does  not  have 
a  significant  effect  on  the  overall  neutral  number  density.  This  does  not  mean,  however,  that  there  is  no 
effect.  Since  the  collision  rate  is  different,  there  might  be  some  difference  that  might  not  be  observable  in 
this  figure. 
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(a)  SCN  At  =  2.5  x  10“7  s 


(b)  SCN  At  =  2.5  x  10“5  s 


Figure  11.  Final  Neutral  Number  Density  for  Coarse  Time  Step 


Figure  12  shows  the  final  number  density  for  the  ions  for  the  coarse  and  fine  time  step  cases.  This  shows 
a  significant  difference  between  the  two  cases.  While  the  main  beam  region  seems  similar,  the  outer  wings  of 
the  plume  are  certainly  not  captured  as  well  in  the  coarse  time  step  case.  Also,  with  the  time  step  so  large, 
the  ions  cannot  make  the  curved  trajectory  to  collide  with  the  back  of  the  ion  source,  which  results  in  the 
difference  between  to  two  cases  in  that  region. 

Figure  13  shows  the  final  electrostatic  potential  for  the  coarse  and  fine  time  step  cases.  This  again  shows 
less  accurate  resolution  of  the  outer  wings  of  the  plume  region.  Also  an  increased  potential  region  behind 
the  plasma  source  exists  since  the  ions  are  not  colliding  with  the  back  of  the  plasma  source  and  becoming 
neutralized. 

Next,  the  local  properties  of  the  plasma  will  be  compared  between  the  three  cases.  The  neutral  sampling 
within  in  the  plume  is  shown  in  Figure  14  for  both  the  velocity  distribution  as  well  as  the  speed  distribution. 
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Figure  12.  Final  Ion  Number  Density  for  Coarse  Time  Step 
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Figure  13.  Final  Electrostatic  Potential  for  Coarse  Time  Step 
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The  most  significant  difference  here  is  that  there  is  a  significantly  larger  band  of  high  velocity  neutrals, 
around  the  20-25  km/s  range,  with  a  noticeable  decrease  in  population  in  the  10  km/s  region  for  the  coarse 
time  step.  The  20  km/s  neutrals  are  results  of  the  charge  exchange  collisions  with  the  high  speed  ions.  The 
10  km/s  neutrals  are  most  likely  from  secondary  collisions  between  the  high  speed  neutral  and  the  other 
particles.  Thus  with  the  large  time  step,  the  high  speed  neutrals  are  able  to  travel  entirely  through  the  high 
density  region  before  another  collision  event  occurs. 
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(a)  SCN  Velocity  At  =  2.5  X  10  7  s 
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(c)  SCN  Speed  At  =  2.5  x  10— 7  s 


(b)  SCN  Velocity  At  =  2.5  x  10“5  s 
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(d)  SCN  Speed  At  =  2.5  x  1CT5  s 


Figure  14.  Neutral  Velocity  and  Speed  Distributions  0.1  m  in  Front  of  Thruster  for  Coarse  Time  Step 


Figure  15  shows  the  ion  sampling  within  the  plume.  Again  a  difference  can  be  seen  between  the  coarse 
and  fine  time  step  cases.  First,  there  is  a  significant  increase  in  the  number  of  ions  in  between  the  two  modes 
in  the  distribution.  This  velocity  range,  around  10  km/s,  corresponds  to  the  decrease  in  the  distribution  of 
neutrals  mentioned  above.  It  appears  that  with  the  coarse  time  step  there  is  not  enough  charge  exchange 
collisions  occurring  between  the  medium  speed  ions  and  neutrals.  Notice  that  these  ions  are  mainly  produced 
by  the  charge  exchange  collisions  from  a  previous  time  step  since  the  plasma  source  is  producing  a  Maxwellian 
distribution  of  ions  with  the  peak  of  the  distribution  at  20  km/s.  Thus  it  appears  that  the  coarse  time  step 
is  moving  the  particles  out  of  the  high  density  region  too  quickly. 

The  neutral  sampling  outside  the  plume  is  shown  in  Figure  16  for  both  the  velocity  distribution  as  well 
as  the  speed  distribution.  While  the  two  distributions  look  similar,  there  is  more  concentrated  collection 
of  particles  in  the  8  km/s  region  for  the  coarse  time  step.  This  is  hard  to  distinguish  from  the  statistical 
scatter  in  the  data,  but  this  clustering  does  occur  over  8  consecutive  velocity  bins,  so  this  may  be  more  than 
a  statistical  artifact. 

Figure  17  shows  the  ion  sampling  outside  the  plume.  In  this  case  there  are  significant  differences  between 
the  coarse  and  fine  time  step  cases.  The  most  notable  is  that  the  coarse  time  step  has  a  bimodal  distribution 
with  the  second  peak  rather  wide  and  centered  around  10  km/s.  This  peak  does  not  appear  on  the  fine  time 
step  and  must  be  a  result  of  the  coarse  time  step.  Again,  the  speed  is  associated  with  secondary  collisions, 
and  if  the  fast  moving  particles,  which  are  the  only  ones  that  create  these  particles,  travel  through  the  high 
density  region  before  another  collision  event  is  performed,  then  there  would  be  left  over  medium  speed  ions. 
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Figure  15. 


Ion  Velocity  and  Speed  Distributions  0.1  m  in  Front  of  Thruster  for  Coarse  Time  Step 
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(b)  SCN  Velocity  At  =  2.5  x  10"5  s 


(d)  SCN  Speed  At  =  2.5  x  10"5  s 

Figure  16.  Neutral  Velocity  and  Speed  Distributions  0.28  m  Above  Thruster  Face  for  Coarse  Time  Step 
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Also,  the  width  of  the  main  peak,  that  is  center  at  4  km/s  is  much  larger  for  the  coarse  time  step  than  for 
the  fine  time  step.  There  is  also  a  corresponding  increase  in  the  speed  distribution  function  value  at  the 
peak.  Again  since  this  peak  is  a  result  of  multiple  collisions  the  coarser  time  step  can  again  attribute  this 
difference. 
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Figure  17.  Ion  Velocity  and  Speed  Distributions  0.28  m  Above  Thruster  Face  for  Coarse  Time  Step 


It  appears  that  there  is  a  significant  coupling  between  the  particle  time  step  and  the  ability  to  capture 
all  of  the  secondary  collisions  that  are  occurring.  Even  though  the  time  step  was  fine  enough  for  the  collision 
modeling  characteristic  time  step,  it  appears  that  the  rapid  particle  density  variation  in  front  of  the  plasma 
source  is  causing  the  collision  characteristic  time  step  calculation  to  be  too  large.  Also,  notice  that  while 
some  of  these  effects  are  noticeable  in  the  gross  perspective  of  the  flow  field,  the  significant  differences  are 
only  apparent  with  the  observation  of  the  velocity  distribution  functions. 

D.  Subcycling  Solution 

Now  that  is  apparent  the  the  coarse  time  step  does  not  capture  a  number  of  significant  flow  features,  an 
analysis  of  the  improvements  associated  with  the  subcycling  algorithm  can  be  performed.  The  subcycling 
algorithm  uses  the  coarse  time  step,  2.5  x  10-5  s,  for  the  slow  particle  time  step  (i.e.,  slow  neutrals)  and  the 
fine  time  step,  2.5  x  1CD7  s,  for  the  fast  particle  time  step  (i.e.,  for  the  fast  ions  and  neutrals).  Table  4  shows 
the  resulting  collision  rates  along  with  the  fine  and  coarse  time  step  collision  rates  for  comparison. 

The  first  thing  to  notice  from  this  figure  is  that  the  subcycling  scheme  results  in  a  decrease  of  computa¬ 
tional  time  from  58.9  hr  to  13.7  hr  compared  to  the  fine  time  step  case.  Unfortunately,  the  collision  rates 
are  nearly  identical  to  the  coarse  time  step  rates  and  are  significantly  lower  than  the  fine  time  step  results 
(with  the  same  exception  of  the  Xe-Xe  elastic  collision  rate). 

While  significant  differences  are  seen  in  the  collision  rates,  Figure  18  shows  the  evolution  of  the  total 
number  of  neutrals  and  ions  as  well  as  the  total  number  of  particles  is  very  similar.  Thus,  the  differences 
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Table  4.  Collision  Rates  for  Subcycling  Case 


Scheme 

Compute 
Time  [hr] 

Time 
Step  [s] 

Total 

Collisions  [#/s] 

Xe-Xe+  Charge 
Exchange  [#/s] 

Xe-Xe 
Elastic  [#/s] 

Xe-Xe + 
Elastic  [#/s] 

SCN 

58.9 

2.5  x  10-7 

5.325  x  107 

5.069  x  107 

2.560  x  105 

2.303  x  106 

SCN 

1.73 

2.5  x  10-5 

4.945  x  107 

4.726  x  107 

2.572  x  105 

1.936  x  106 

SCY 

17.3 

2.5  x  10-5 

4.966  x  107 

4.745  x  107 

2.535  x  105 

1.963  x  106 

between  the  subcycling  case  and  the  fine  time  step  case  again  must  be  for  only  a  small,  but  significant, 
fraction  of  the  total  particles. 
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Figure  18.  Global  Particle  Counts  for  Subcycling  Case 


Figure  19  shows  the  final  number  density  for  the  neutrals  for  the  subcycling  case  as  well  as  the  fine  and 
coarse  time  step  cases.  This  case  again  shows  that  there  is  little  difference  between  the  neutral  number  den¬ 
sities  for  the  subcycling  case.  Thus,  the  overall  neutral  distribution  is  fairly  insensitive  to  the  computational 
time  step.  However,  as  was  the  case  for  the  coarse  time  step,  it  is  expected  that  there  will  be  differences 
that  are  not  observable  in  this  perspective. 


Figure  19.  Final  Neutral  Number  Density  for  Subcycling  Case 


Figure  20  shows  the  final  number  density  for  the  ions  for  the  subcycling  case  as  well  as  the  fine  and 
coarse  time  step  cases.  Unlike  the  coarse  time  step  case,  the  ion  number  density  is  very  close  to  the  fine  time 
step  case.  The  outer  wings  of  the  plume  are  captured,  and  the  ion  neutralization  at  the  back  of  the  plasma 
source  is  also  captured.  Therefore,  the  subcycling  is  drastically  improving  the  capabilities  of  capturing  the 
ion  distribution. 

Figure  21  shows  the  final  electrostatic  potential  for  the  subcycling  case  as  well  as  the  fine  and  coarse  time 
step  cases.  This  again  shows  that  the  subcycling  case  and  the  fine  time  step  cases  are  quite  similar.  This  is 
to  be  expected  since  the  electrostatic  potential  is  directly  related  to  the  ion  distribution.  Even  the  potential 
drop  behind  the  thruster  seen  in  the  fine  time  step  case  is  captured  in  the  subcycling  case. 

While  the  collision  rates  are  different,  the  subcycling  case  has  so  far  improved  the  ion  number  density 
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Figure  20.  Final  Ion  Number  Density  for  Subcycling  Case 


Figure  21.  Final  Electrostatic  Potential  for  Subcycling  Case 
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distribution  compared  to  the  coarse  time  step  case  and  has  shown  no  difference  in  the  neutral  number  density 
distribution. 

Next,  the  local  properties  of  the  plasma  will  be  compared  between  the  subcycling  case  and  the  fine  and 
coarse.  The  neutral  sampling  within  in  the  plume  is  shown  in  Figure  14  for  both  the  velocity  distribution 
as  well  as  the  speed  distribution.  Unfortunately,  the  subcycling  cases  looks  much  more  similar  to  the  coarse 
time  step  case  and  has  significant  differences  with  the  fine  time  step  case.  The  same  arguments  about  the 
coarse  time  step  differences  also  seem  to  apply  here.  While  the  neutrals  are  propagating  at  the  fine  time 
step,  there  is  still  no  mechanism  to  get  these  neutrals  to  participate  in  collision  events  while  they  reside  in 
the  high  density  regions. 


(a)  SCN  Velocity  At  =  2.5  x  10  '  s 


speed  (m/s) 

(d)  SCN  Speed  At  =  2.5  X  10“7  s 

Figure  22. 


(b)  SCN  Velocity  At  =  2.5  x  1CT5  s 


speed  (m/s) 

(e)  SCN  Speed  At  =  2.5  x  10“5  s 


(c)  SCY  Velocity  At  =  2.5  X  10  5  s 
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(f)  SCY  Speed  At  =  2.5  x  1CT5  s 


Neutral  Velocity  and  Speed  Distributions  0.1  m  in  Front  of  Thruster  for  Subcycling  Case 


Figure  23  shows  the  ion  sampling  within  the  plume.  This  case  also  shows  significant  differences  between 
the  subcycling  and  fine  time  steps  and  is  very  similar  to  the  coarse  time  step.  It  appears  that  the  secondary 
collision  discussion  from  above  explains  these  differences. 

The  neutral  sampling  outside  the  plume  is  shown  in  Figure  24  for  both  the  velocity  distribution  as  well 
as  the  speed  distribution.  While  it  is  not  certain  that  the  8  km/s  region  is  the  coarse  time  step  is  caused  by 
statistical  scatter,  it  is  worth  noting  that  the  subcycling  case  does  not  demonstrate  this  feature.  Otherwise, 
the  subcycling  case  looks  very  similar  to  the  fine  time  step  case. 

Figure  25  shows  the  ion  sampling  outside  the  plume.  This  case  shows  drastic  improvements  from  the 
coarse  time  step.  The  subcycling  case  does  not  have  the  secondary  peak  in  the  10  km/s  range  and  has 
a  similarly  narrow  speed  range  around  the  most  probable  speed.  It  is  apparent  that  the  subcycling  does 
significantly  improve  the  particle  modeling  outside  the  high  density  plume  region.  This  is  most  likely  due 
to  the  fact  that  the  trajectory  of  the  high  speed  ions  is  significantly  improved  with  the  fine  time  step,  and 
thus  these  high  speed  ions  are  more  effected  by  the  electrostatic  potential  field. 

IV.  Conclusions 

The  acceleration  subcycling  acceleration  scheme  within  the  AQUILA  plasma  modeling  module  of  COL¬ 
ISEUM  was  investigated  to  determine  how  effective  it  is  in  capturing  the  local  plasma  properties.  First, 
the  simulation  was  demonstrated  to  be  capable  of  converging  to  a  solution  for  a  sufficiently  fine  time  step. 
This  solution  was  the  used  to  compare  the  performance  of  the  simulation  at  a  much  coarser  time  step.  This 
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(a)  SCN  Velocity  At  =  2.5  x  10  '  s 


(b)  SCN  Velocity  At  =  2.5  x  10  5  s 


(c)  SCY  Velocity  At  =  2.5  x  10  5  s 


Figure  23.  Ion  Velocity  and  Speed  Distributions  0.1  m  in  Front  of  Thruster  for  Subcycling  Case 
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Figure  24. 


Neutral  Velocity  and  Speed  Distributions  0.28  m  Above  Thruster  Face  for  Subcycling  Case 
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Figure  25.  Ion  Velocity  and  Speed  Distributions  0.28  m  Above  Thruster  Face  for  Subcycling  Case 


showed  several  deficiencies  in  the  coarse  time  step  solution.  These  were  mainly  focused  on  the  fact  that  the 
high  speed  particles  are  leaving  the  high  density  region  where  multiple  collisions  are  expected  to  occur  after 
one  or  two  time  steps.  In  addition  the  coarse  time  step  resulted  in  the  ion  trajectories  being  significantly  off 
which  resulted  in  large  differences  in  the  ion  density  distribution  compared  to  the  fine  time  step  case. 

The  subcycling  case  improved  the  modeling  of  the  number  densities  of  the  ions  and  neutrals  compared  to 
the  coarse  time  step  simulation  with  a  compute  time  speedup  of  a  factor  of  3.4.  It  also  significantly  improved 
the  modeling  of  the  lower  density  region  outside  the  main  plume.  In  this  region  the  velocity  distribution 
functions  for  the  neutrals  and  ions  were  both  very  similar  to  the  fine  time  step  case.  This  shows  that  by 
improving  the  modeling  of  the  electrostatic  forces,  via  finer  time  steps  to  propagate  the  ions  and  the  force 
calculations,  does  improve  the  modeling  of  the  plasma.  For  the  higher  density  region  of  the  main  beam 
of  the  plasma,  the  subcycling  showed  no  improvement  compared  to  the  coarse  time  step  solution.  While 
this  appears  to  be  rooted  in  the  fact  that  a  high  speed  particle  can  travel  through  the  entire  high  density 
plume  region  in  one  time  step  and  thus  not  participate  in  the  additional  collisions  that  appear  to  produce 
the  different  features  of  the  fine  time  step  solution.  An  alternative  explanation  to  this  behavior  could  be 
that  the  collision  model  is  not  providing  consistent  results  at  the  different  time  scales.  Suppose  the  collision 
model  was  selecting  too  many  collisions  to  occur  at  the  finer  time  step  because  the  selection  criteria  was 
slightly  off.  This  would  result  in  a  non-physical  increase  in  additional  collisions  that  the  high  speed  particles 
would  be  participating  in.  This  would  mean  that  the  fine  time  step  solution  was  not  correct. 

To  further  identify  the  root  cause  of  this  problem,  it  is  apparent  that  the  collision  model  with  AQUILA 
needs  to  be  further  validated  in  order  to  assure  that  it  is  not  producing  a  surplus  of  collision  events  as  the 
time  step  is  decreased.  While  the  collision  model  has  been  reported  to  have  been  successful  in  previous 
studies,2,3,8, 14  a  systematic  validation  of  the  collision  model  source  code  has  not  been  found  by  the  authors. 

Another  suggested  improvement  to  address  the  remaining  differences  between  the  subcycling  and  the 
fine  time  step  cases  is  to  perform  the  collision  modeling  between  the  high  speed  particles  and  the  rest  of 
the  particles  within  the  subcycle  loop.  While  this  would  significantly  decrease  the  computational  efficiency 
associated  with  this  scheme,  it  should  address  the  problems  in  the  high  density  regions.  Performing  this 
modification  would  be  a  non-trivial  task  as  the  collision  code  would  need  to  be  modified  to  only  perform 
certain  collision  pairing  at  certain  times. 
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